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Abstract 



We study properties of two-leg Heisenberg spin ladders in a mean-field ap- 
proximation using a variety of angular-momentum coupled bases. The mean- 
field theory proposed by Gopalan, Rice, and Sigrist, which uses a rung ba- 
sis, assumes that the mean-field ground state consists of a condensate of 
spin-singlets along the rungs of the ladder. We generalize this approach to 
larger angular-momentum coupled bases which incorporate — by their mere 
definition — a substantial fraction of the important short-range structure of 
these materials. In these bases the mean-field ground-state remains a conden- 
sate of spin singlet — but now with each involving a larger fraction of the spins 
in the ladder. As expected, the "purity" of the ground-state, as judged by the 
condensate fraction, increases with the size of the elementary block defining 
the basis. Moreover, the coupling to quasiparticle excitations becomes weaker 
as the size of the elementary block increases. Thus, the weak-coupling limit 
of the theory becomes an accurate representation of the underlying mean- 
field dynamics. We illustrate the method by computing static and dynamic 
properties of two- leg ladders in the various angular- momentum coupled bases. 



I. INTRODUCTION 



Spin ladders have received much recent attention due to their likely relevance to high- 
ly superconductors. The full arsenal of theoretical many-body physics methods has been 
brought to bear on the problem of spin ladders, including direct diagonalization (DD) [|TJ, 
Lanczos techniques ||, quantum Monte Carlo (QMC) density matrix renormalization 
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group (DMRG) 0, high-order perturbation theory augmented by re-summation techniques 
HI, mappings to effective field theories [|||| as well as quasi-analytic methods based on mean- 
field approaches . In the present work, we focus on the last of these techniques and extend 
the mean- field theory (MFT) of Gopalan, Rice and Sigrist J7|. This extension includes a 
derivation of generalized mean-field equations appropriate to any angular-momentum cou- 
pled basis. In particular, we examine MFT's utilizing bases whose elements are eigenstates 
of 2x1 ( "rung" basis) , 2x2 ("plaquette" basis) and 2x4 systems. The original work of Ref. J7| 
employed a rung basis. Indeed, the strong coupling limit for which coupling within rungs is 
much greater than coupling between rungs — in which case the ladder ground state can be 
well-approximated by a configuration in which every rung is coupled to a spin singlet — was 
the inspiration for their MFT in which rung singlets "condense" . We also do MFT calcula- 
tions in the rung basis but obtain rather different results from those of Ref. j7|. The likely 
origin of these differences is addressed below. We find our rung-basis MFT results in much 
better agreement with known properties of 2-leg ladder — such as the ground state energy 
per site, the spin gap, the spin correlation length ||, the one-magnon dispersion relation 
||, and the strength of the S(n, ir) spin response || — than those of Ref. |7|]. Our results 
suggest that rung-basis MFT is a relatively simple but still robust means of accounting for 
the physics of spin ladders. 

We also study plaquette- and 2x4-basis MFTs. In a previous work ||10|| , we showed that 
the plaquette basis offers some significant advantages over simpler bases (including the rung 
basis) when describing the physically relevant JIT | case of isotropic spin ladders where intra- 
rung and inter-rung couplings are equal. We find here that the plaquette basis yields an 
MFT which is superior to that of the rung basis. While the two sets of calculations give very 
similar mean-field values for many ladder observables, we conclude that only the plaquette 
basis results are reliable. This determination is based on the fact that our estimates of 
contributions beyond mean-field theory are quite large for the rung basis but much smaller 
for plaquettes. As our general approach holds for any angular-momentum coupled basis, we 
also study the MFT utilizing a 2 x 4 basis. Here, MFT results are comparable to those for 
the rung and plaquette MFTs. We trace the lack of significant improvement upon going to 
this rather complicated basis to truncations which become increasingly severe as the basis is 
built from the eigenstates of larger and larger systems. In this context, the plaquette basis 
MFT seems the best compromise between building enough physics into the basis that MFT 
can account for the bulk of the physics and, at the same time, avoiding severe truncations. 

The paper is organized as follows: In Section |IJ, we present our angular-momentum 
coupled formalism for spin ladders in the context of the rung basis. We also present our 
rung-basis MFT and identify differences with Ref. M. In addition, we derive weak-coupling 
expansions for the energy per site and spin gap which are useful in understanding the 
physics embodied in the MFT equations (and which additionally reproduce the full MFT 
results with surprising accuracy, especially for the plaquette basis). We next indicate how to 
make a perturbative estimate of corrections to the energy per site and spin gap arising from 
terms in the Hamiltonian which are omitted at the mean-field level. Finally we extend our 
MFT to plaquette and 2x4 bases and indicate where and with what justification truncations 
are made. Numerical results and discussion of them appear in Section [TTI|, while Section [TV] 
contains our conclusions. 
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II. FORMALISM 



A. Two-leg Ladder Hamiltonian 

The Heisenberg Hamiltonian for a two-leg ladder with nearest-neighbor antiferromagnetic 
coupling is given by 

H=J\\Y,Si-S j + J±Y, S i- S J- C 1 ) 

Here J± and J\\ are the strength of the antiferromagnetic coupling along the rungs and 
chains of the ladder, respectively 0. In the strong-coupling limit, A = J\\/ J± 0, the 
interaction between rungs is very weak and the above separation of the Hamiltonian is, 
indeed, a natural one. In this limit the Hamiltonian can be separated into a large component 
H , which contains the strong intra-rung interactions, and a weak "two-body" potential V 
between rungs. That is, H = H + XV, where 

N 

tf = ESi(r)-S 2 (r), (2) 

r=l 

N 

V = J2 [Si(r) • Sx(r + 1) + S 2 (r) • S 2 (r + 1)] . (3) 



r=l 



Note that we have set J± = 1 and that r labels a specific rung in the A^-rung ladder. Further, 
Si(V)[S2(r)] denotes the spin operator for the electron on the first [second] site of the r-th 
rung. The rung basis is particularly suitable in the strong-coupling regime, as H is diagonal 
in this basis. The eigenvalues and eigenvectors of H are of the form: 

\<f>t m ) = I (s lS2 ) im) ■ q = U(£ + 1) - 3/4 . (4) 

Here the two spins along a definite rung in the ladder are coupled to a total angular mo- 
mentum £, which can be zero or one, and projection m. Rung singlet and triplet creation 
operators are defined accordingly 



c\(r) <g> c 2 (r) 



em 



(5) 



where is an electron creation operator. Note that the rung singlet/triplet operators satisfy 
the usual commutation relations for bosons: 

[Ai m (r), A], m ,(r')] = 8u'$mm'$rr' ■ (6) 

We now proceed to rewrite the above Hamiltonian in a second-quantized form. For the 
unperturbed part of the Hamiltonian, namely, the part that is diagonal in the rung basis, 
we obtain 

N 1 N 

Ho = EE^[4(0®^)] 00 = E [e sHr)s(r)-e 1 t\r)-t(v)] , (7) 

r=l£=0 ' r=l 
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where we have defined t = ^/2£ + 1 and 



S t(r)=4 (r); 4(r) = 4» ■ (8) 
Note that we have also introduced the "tilde" operators through the definition: 

tm ( 1) t— m ■ (9) 



As we have indicated explicitly in Eq. (|7]) it is the operator t — not t — that transforms as 
a vector under rotations; the m-dependent phase factor in the above equation is necessary 
to enforce this condition. In analogy to the one-body term, the two-body potential energy 
[Eq. (|3])] can be written in a rotationally-invariant form 



N 



v = J2 £ 5(^11^(1,2)11^') 



l,(r)Al,(r + l)).(A £l (r)A h (r + l)) : 



oo 



(10) 



Matrix elements of the potential have been computed in Ref. || and have been listed — 
along with the corresponding matrix elements in other bases — in Table Q. It is important to 
mention that due to the symmetry of the potential only two of these matrix elements are 
independent. We have chosen them to be: 



a= (011i|V(l,2)||101) and = <110||V(1, 2)||110) . 



(11) 



If we now assume periodic boundary conditions, then the total linear momentum be- 
comes also a good quantum number. Thus, it is convenient to re-write the Hamiltonian in 
momentum space. To this end, we start by performing a Fourier transformation of the spin 
singlet and triplet operators 



£ e%kl Aem(^) 



k 



2nu 



u = l,2 N 



Moreover, we decomposed the Hamiltonian into three components: 

H = Hq + XH\ + XH2 ■ 



(12) 



(13) 



The first term in this decomposition is simply the unperturbed Hamiltonian of Eq. (|7|) 
written in momentum space 



H = ^2[^(k)s(k)-e 1 t\k)-t(k) 



(14) 



The second term represents that part of the two-body potential that is bilinear in both the 
spin-singlet and spin-triplet operators. It is given by 

# 1 = ^ £ a cos(g) \t\k[) ■ t j! (k' 2 )s(k 2 )s(k 1 ) - t\k[) ■ t(A; 2 )s t (^)s(A; 1 ) + h.c.j . (15) 



k\k2q 



Note that we have introduced the two "final" momenta by defining k[ = k±— q and k' 2 = k 2 +q. 
Finally, H 2 is the remaining — quartic in the spin-triplet operator — part of the potential: 
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j=0 k±k 2 q 



00 



(16) 



where Cj = (— l) J+1 j (1— 3<5j-o)/2. Note that there is no component of the potential having an 
odd number of spin-triplet operators because of the selection rule l\ + £ 2 + £[ + £' 2 = even || . 

We now conclude this section by introducing the "grand-canonical" Hamiltonian K 
through the definition 

K = H — jjN . (17) 
Here \x is the chemical potential and N is the number operator given by 

N = Y, [s\k)s(k) + £ tl(k)t m (k)} = £ P (k)s(k) - t\k) ■ t(k)} . (18) 



As long as no approximations are being made, the description of the interacting system 
in term of either H or K is equally acceptable. The virtue of K, however, is in that it 
allows a consistent treatment of the system even if the "particle" number is not conserved 
(i.e., if [H,N] ^ 0). Since the mean- field approximation of the next chapter is based on 
a Bogoliubov transformation, using the grand-canonical Hamiltonian K provides a definite 
advantage over H. 



B. Mean-Field Approximation 

In this section we implement the mean-field approximation of Gopalan, Rice, and 
Sigrist H . The central dynamical assumption of their model is that the spin-singlet bosons 
with zero linear momentum form a condensate. Hence, it becomes natural to replace the 
spin-singlet operators by "c-numbers", i.e., 

s{k),s\k) — >VNs8 kt0 . (19) 



If this condition is now substituted on the above expressions for H and Hi [Eqs. ( |I4| , |T5| ), 
respectively] we obtain 

H = Ns 2 e -e 1 ^(k)-t(k) , (20a) 

k 

H 1 = s 2 Y J u cos{k) [t\k) ■ t^-ife) - tf(k) ■ t{k) + h.c] . (20b) 

k 

These expressions are now bilinear in the spin-triplet operators and, thus, amenable for 
diagonalization via a Bogoliubov transformation. Note that H 2 [Eq. fll6f) 1 remains quartic 
in the spin-triplet operators. Hence, for the moment we will ignore this term in the de- 
termination of the mean-field ground state and will return to it later to compute its effect 
perturbatively. In this way, the mean-field Hamiltonian can be written in the following 
compact form 

K MF (s 2 , /x) = A(e - ^)s 2 - |jV( Cl - M) + \ E T K k ) V(k) T m (k) . (21) 

km 
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where 



t m {+k) 



and Q(k) 



A(k) -2A(k)\ 
-2A(k) A(k) J 



.4.H0. 

Note that the two independent elements of the 2x2 matrix fi(fc) are given by 
A(k) =£]_ — // + 2Xas 2 cos(fc) and A(A;) = Xas 2 cos(fc) . 



(22) 



(23) 



The above mean-field Hamiltonian can now be diagonalized by means of a standard Bogoli- 
ubov transformation. We obtain 



K MF (s\ /,) = N(e -fi)s 2 - ^N(e l - y) + | £ u(k) + £ u(k)jl {k)-y m (k) 

k km 



where 



uj(k) = \j A(k) 2 - 4A 2 (fc) = (ei - p)^Jl + dcos{k) ; (d = 4Aas 2 /(ei - p) 
Moreover, the new "quasiparticle" creation operators are defined by 

>yl(k) = u(k)tl(k) - v(k)t m (-k) , 



(24) 



(25) 



(26) 



where the coefficients generating the canonical transformation — satisfying [u 2 (k) — v 2 (k)} 
1 — are given by 



u(k) 



A(k)+Lj(k) 

\\ 2u(k) 



and v(k) = sgn[A(fc)], 



A(fc) - u(k) 
\\ 2u{k) 



(27) 



The mean-field ground state — the quasiparticle vacuum — is now determined by the condition 

7m (fc)|$ > = , for all k and m . (28) 

As we have described in the appendix, we have carried out the Bogoliubov transformation by 
diagonalizing an appropriate 2x2 random-phase- approximation (RPA) matrix. Although 
not explicitly treated in the appendix, the method has the virtue of being easily generalizable 
to matrices of arbitrary dimension. We should mention that for the particular case of the 
rung basis (i.e., e = — 3/4, e\ = 1/4, and a = 1/2) our results are identical to those obtained 
in Ref. — with one exception. The two factors of 3/2 appearing in Eq. ( p4|) appear as 
factors of 1/2 in Ref. 0. In our derivation these factors arise from the zero-point motion of 
the three triplet "oscillators", one for every magnetic component of t. We will examine the 
quantitative impact of this difference in Sec. pi . 

Ground-state observables can now be obtained by taking appropriate derivatives of the 
thermodynamic potential 

K (s 2 ,fi) = ($ |#mf(s 2 ,//)|$o> = N(e - fi)s 2 - -A^ - y) + ^#) . (29) 

~ k 



For example, the number of rungs can be recovered by using the thermodynamic relation 
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N = — 



'dK (s 2 ,/j,)' 



(30) 



or equivalently, from the following transcendental equation 



2 5 3 r 

s 2 h - 

2 2 



1(d) - dl'(d) = ; 1(d) = 



27r dk 
o 27 



1 + dcos(k) 



(31) 



This equation serves to write d — and, thus, the chemical potential \x — in terms of s 2 . Note 
that the function d(s 2 ) is "universal"; it is independent of the value of the ratio of Jn to J 



_L 



and of the choice of basis. Finally, the expectation value of the singlet condensate may be 
determined from minimizing the ground-state energy with respect to s 2 



'dE n (s 2 ) 



ds 



o2 



0; E (s 2 ) = K (s 2 ,v(s 2 ))+v(s 2 )N 



(32) 



Having determined the values of the chemical potential \x and the singlet condensate s 2 , static 
observables — such as the energy-per-site and the one-magnon dispersion relation ou(k) — can 
be easily computed. 



C. Dynamic Spin Response 

Perhaps more exciting is the possibility of computing the dynamic spin response of the 
system. This kind of computation has proven to be particularly challenging in certain 
approaches — such as DMRG — as transition matrix elements, in addition to the excitation 
spectrum, must be computed. Yet the dynamic response constitutes one of the most fun- 
damental properties of the system and one that may soon be extracted from new data on 
inelastic neutron scattering experiments from single crystals. 

For two leg-ladders one can define two independent spin responses, So(q, uj) and S n (q, cu), 
according to the relative phase between the two spin operators along the rungs of the ladder. 
That is, 



S ,«(q,u) = £ktf„|<V(5)|*o) S(u-u n ) , 



(33) 



where \l/o is the exact ground-state wavefunction of the system, \l/ n an excited state with 
excitation energy u> n , and the transition operator is given by 



^(g) = E e ^ r f s iW±S 2 (r 



(34) 



Although similar in form, these two operators probe very different aspects of the dynamics 
of the system. This can be seen most easily by expressing both transition operators in a 
second-quantized form: 



Si(r) + S 2 (r)l =v / 2[t t (r)®t 



lm 



Si(r)-S 2 (r) 



t^(r)s(r) - s f (r)t m (r) 



MF 



tl(r) - t m (r) 



(35a) 
(35b) 
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These equations imply that S n probes the low-energy part of the spectrum, as it couples 
the ground-state of the system to the one-magnon band. In contrast, Sq — a two-magnon 
operator — probes the high-energy part of the response. Thus, the S n response becomes 
instrumental in elucidating the important low-energy/low-temperature behavior of the sys- 
tem. 

Matrix elements of the above operators can be readily evaluated by expressing the spin- 
triplet operators in terms of the quasiparticle operators of Eq. (BBT) . We obtain 



N 
1 



S (q,u) 



2 * dk 
o 2^ 



u(k)v(k + q) — u(k + q)v(k) 8{uj—uj{k)—uj{k + q)\ ; 



S v (q,u) = s 2 u(q) -v(q) 6(u-u(q) 



s 2 



1 + dcos(q) 



5(uj—uj 



(36a) 
(36b) 



Note that, given q and u, the delta function in Eq. ( |36a|) constrains the integral over k 
with the result that not only is the Sq response pushed to high-energy but, in addition, 
it is strongly fragmented. This is in contrast to the SV response that is predicted to be 
concentrated in a single excitation at uj—u(q). At least for the S n (q = 7T,uj) response, this is 
in good agreement with an earlier calculation |J on a 2 x 16 ladder that predicts almost 95 
percent of the strength to be concentrated at an excitation energy equal to the singlet-triplet 
gap. Finally, having computed the two dynamic response functions the spin-spin correlation 
function per site can be easily extracted: 



($o|5 2 (0)5 2 (r)|$ ) = I [** p. e *\S (q) + S v {q)]/N 

8 Jo 2tt l j 

where the static structure factors have been defined by 



Jo 



(37) 



(38) 



Note that in Eq. (|37|) both spins are located on the same chain of the ladder at a distance of 
r rungs. As we will see shortly the long-range behavior of the spin-spin correlation function 
is dominated by the static structure factor S 7T (q), namely, 



($ o \S z (0)Sz(r)\®o) 



i7r dq 
— < 
2tt 



■iqr 



s 2 



i7r dq 



2tt 



dcos(g) 



(39) 



We can see that, in the spirit of the stationary-phase approximation, the principal contri- 
bution to the integral at large r occurs where the rapid oscillation of the exponential factor 
is canceled by rapid changes in [1 + <icos(g)] -1 / 2 . This will happen near any of the singu- 
larities of the latter factor. Indeed, upon using the steepest-descent technique we obtain 
the expected exponential behavior of the spin-spin correlation function and our analytic 
approximation for the correlation length: 



(%\S z (0)S z (r)\%) 



~ e 



£ = 1/ln d/{\ - VI - d 2 ) 



(40) 
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D. Weak-Coupling Expansions 



To gain further insight into the mean-field theory of Ref. @ we derive in this section 
weak-coupling (A < 1) expansions for various observables. We start this procedure with 
an expansion of Eq. (j31~l), which determines the parameter d — or equivalently the chemical 
potential — in terms of the parameter P = l — s 2 . That is, 



d 



45 / d 2 



16 



16 



-t 2 



(41) 



This equation can now be substituted into the ground-state energy to yield an expression 
that is correct to fourth-order in the small parameter t = O(X): 



E Q {t) 



St + A P - -t 3 + 0{&) ; where 5 = VdX\ 
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(42) 



By minimizing the above expression for the energy with respect to i, ground-state — and 
even some excited-state — observables can be determined. Indeed, from the minimization 
procedure the singlet condensate is readily obtained 



s 2 



3 /AaY 
~2 I ^ J 



(43) 



Moreover, having determined the behavior of the singlet condensate in terms of A, one can 
generate — through Eqs. ([25]) and (4"i"|) — the corresponding behavior of the chemical potential. 
This is all that is required to compute the mean-field observables. For example, we obtain 
the following perturbative expansions for the ground-state energy per rung and for the 
singlet-triplet gap, respectively: 



En 



A 



3(Aa) 



2 A 
A - 2A|a| 



3 / Xa 
- ( X\a\ 



An 



+ 



C(A 4 ) 
1 / Ala 



An 



0{X 3 



(44a) 
(44b) 



Note that the small parameter in the above expansions is given by Aa/A . 



E. Perturbative Contribution From H2 

An important assumption of the mean-field approximation of Ref. J7| is that H2 — the 
term of the Hamiltonian quartic in the spin-triplet operators [Eq. (|16"|)1 — changes the mean- 
field results only slightly. To test the validity of this assumption we now estimate the effect 
of H 2 on the ground-state energy and on the one-magnon dispersion relation using first-order 
perturbation theory. 

The first-order correction to the ground-state energy due to H 2 is given by 

A£ = A($o|# 2 |$o> • (45) 



9 



To evaluate this correction to the energy the four spin-triplet operators in H2 [Eq. ([TED ] 
are expanded in terms of the quasiparticle operators — J m (k) and ^(k) — and the resulting 
matrix elements are evaluated using Wick's theorem. We obtain, 



AE /N = 3A/3 



27r dk \ 2 

— u(k) cos(k)v(k) 
2n ' 



2n dk \2 

— v(k) cos(k)v(k) 
2ir 



(46) 



Note that only those combinations of quasiparticle operators that conserved quasiparticle 
number need to be considered. The corresponding perturbative contribution to the one- 
magnon band is slightly more complicated to evaluate. Still, it can be efficiently computed 
from the following double-commutator formula: 

Au(k) = \[(<£> \lm(k)H 2l l(k)\$o) ~ ($q|# 2 |$o>] = \(*o\[y m {k), [H 2 ,yUQ]]\*o) • (47) 

In this way the resulting expression for the first-order correction to the dispersion relation 
becomes 



Au(k) = 2\p\2u(k)v(k) cos(k) 



ln dq 
2^ 



u(q) cos(q)v(q) 



2tt 



dq 



u (k)+v (k) cos(k) / — v(q) cos(q)v(q) 

J JO 27T 



(48) 



F. Mean-Field Theory with Improved Bases 

The mean-field theory of the previous section is fully characterized by three parameters: 
the spin-singlet energy 6q, the spin-triplet energy ei, and the singlet-triplet matrix element 
a (if corrections to the mean-field observables need to be computed, an additional triplet- 
triplet matrix element, /3, is necessary). This suggests that the generalization of the mean- 
field method of Gopalan, Rice, and Sigrist to larger angular-momentum coupled bases is 
straightforward. Moreover, the advantage of the larger bases is that they incorporate part 
of the short-range structure of the system — and, thus, an important part of the physics— 
exactly. 

One such basis is the "plaquette" basis introduced by us in Ref. fLOf . The plaquette basis 



is generated from coupling the four spins in a 2 x 2 lattice. Specifically, the two diagonal 
pairs of spins are coupled to well-defined total angular momenta, £14 and £23, which can equal 
zero or one. These two "link" angular momenta are in turn coupled to a total plaquette 
angular momentum j and projection m. The physically appealing feature of the plaquette 
basis is that the 2x2 ladder Hamiltonian with isotropic (Jn = Jj_ = 1) coupling is diagonal 
in the plaquette basis. That is 

#(2x2) t u l&,jm) = + 1) - e u (£ u + 1) - £ 23 {£ 23 + 1)1 kd^jm) . (49) 



The lowest spin-singlet (j — 0) eigenstate arise from coupling the two link angular momenta 
to their maximum value of £14 = £23 — 1. The corresponding value for the energy is given 
by e = —2. Next comes the three-fold degenerate spin-triplet (£ u = £23 = j = 1) eigenstate 
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with an energy of e\ = —1. Note that these four states — out of a total of sixteen — define 
the low-energy basis the we have used earlier in a study of the dynamic spin response ||. 
It is precisely these four low-energy states that we retain here to generalize the mean-field 
theory of Ref. [7| to the plaquette basis. The sole remaining parameter needed to define the 
mean-field theory is the singlet-triplet matrix element a. In the plaquette basis it is given 
by a = —1/3. It is interesting to note that the ratio of a to the unperturbed gap Ao — which 
is a measure of the "goodness" of the basis [see Eq. ([P])] — improves as one goes from the 
rung to the plaquette basis: |a/A | goes from 1/2 to 1/3. 

We can also compute the two dynamic spin responses in the plaquette basis. Indeed, 
the So and S n operators of Eq. (|33|) can be re-expressed in terms of the four spin operators 
along a 2 x 2 plaquette. We obtain 



S (q) =e- ig / 2 ^2e l2qp 



S w (q) = e -n/2J2- i2qp 



e 



cos(g/2)(Si(p) + S 2 (p) + S 3 (p) + S 4 (p) 
isin(g/2)(S 1 (p) + S 2 (p) - S 3 (p) - S 4 (p) 
cos(g/2)(s 1 (p) - Sa(p) + S 3 (p) - S 4 (p) 
-isin(g/2)fS 1 (p) - S 2 (p) - S 3 (p) + S 4 (p) 



(50a) 



(50b) 



Note that Si(p)[S 3 (p)] and S 2 (p)[S 4 (p)] denote the first and second spin operators along 
the first [second] rung of the p-th plaquette. Moreover, the sum over p now runs over all 
plaquettes in the ladder and the Fourier sum includes 2q — not q — as the variable conjugate 
to p. In spite of their apparent complexity the above expressions are identical to the ones 
given in Eq. (|4]). The virtue of writing them in such a form becomes evident as one takes 
their mean- field limit. Indeed, in this limit the dynamic spin responses take the following 
simple form: 

— S (q,u) = ; (51a) 



N 
1 



^- r S^q,u) = -s 2 sin2 ^/ 2 ) = sL-u(2q)) . (51b) 
l + dcos(2g) v 1 



Perhaps more than anything else these expressions capture the essence of the mean-field 
approximation in the plaquette basis. Clearly, the S response is not zero [see. Eq. Q36a|) ]. 
Yet, it is small and its strength is fragmented and located at very high excitation energy. 
The "low-energy" plaquette basis — the truncated basis obtained by retaining only the lowest 
spin-singlet and spin-triplet states — can not account for the physics at high excitation energy. 
In contrast the S n response, which couples the ground-state to the one-magnon band, probes 
the low-excitation part of the spectrum and should be well described using the low-energy 
plaquette basis. 

One could continue to systematically improve the basis. For example, the spectrum of 
isotropic 2x4 and 2x8 ladders have been previously calculated in Ref. ||. In particular, the 
lowest singlet and triplet energies of the 2x4 ladder (with open boundary conditions) are 
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given by eo = —4.29307 and e\ = —3.52286, respectively. Moreover, from the 2x8 calculation 
the singlet-triplet matrix element can also be extracted; it is given by: a = — 0.20805^—1/5. 
Thus, the mean-field theory can be further extended to the 2x4 basis. Note that for this 
basis |a/A | ~ 1/4. Although one has obtained a better expansion parameter in this basis 
relative to the plaquette basis, the improvement has been somewhat reduced by the fact that 
as the singlet-triplet matrix element a goes down so does the unperturbed gap. Moreover, 
the truncation errors increase considerably; while one retains the same four states in both 
bases, the dimension of the vector space goes from 16 in the plaquette basis to 256 in the 2x4 
basis. This will become an important issue in our final analysis. Values for the mean-field 
parameters in the various bases are tabulated in Table |. 



III. RESULTS 

In this section we report mean-field results for the energy-per-site, the singlet-triplet gap, 
the one-magnon dispersion relation, the spin-spin correlation function, and the dynamic spin 
response using various angular-momentum coupled bases. Because of the nature of these 
bases, any comparison among them will only be carried out at the isotropic (A = 1) point. 



A. Rung Basis 

We start with a presentation of our results in the rung basis. In Table [TT] we have listed 
the ground-state energy-per-site for three values of A. The second column contains the 
mean-field numbers reported in Ref. J7J, while the third column contains the corresponding 
values as obtained by us. The difference between these values originates solely from the 
treatment of the zero-point motion; recall that a factor of 1/2, rather than 3/2, was used 
in Ref. J7j. While the impact of this error on the spin gap has been discussed recently in 
Ref. [12 1, we feel compelled to address it here as well, as the mean-field theory of Gopalan, 



Rice, and Sigrist seems to work much better for most observables — not only for the spin 
gap — than their original results might suggest. Indeed, even for the most unfavorable case 
of A = l, the corrected mean-field results are within 6 percent — not 18 percent — of the exact 
answer. Fig. [I] presents these results in graphical form and serves as a further testimony to 
this fact. 

For the spin gap the improvement is even more dramatic. In Table H and Fig. ^| we display 
results for the singlet-triplet gap (in units of J±) as a function of A. The filled circles in 
the figure (labeled "exact" ) were obtained from a very high-order perturbative expansion Q 
supplemented by a Pade analysis. This plot shows conclusively that the mean-field theory of 
Ref. captures an important part of the physics governing the ladder materials. While the 
perfect agreement at the isotropic point seems a bit fortuitous, it is clear that the mean-field 
theory provides — in addition to an intuitive physical picture — a very robust starting point 
for more sophisticated calculations. 

In Fig. |3] we show our mean-field results for the spin-triplet dispersion band (in units of 
J±) relative to the value at the band minimum A = uu(jr). The filled circles for A = 0.1 and 
0.5 are the result of a perturbative calculation H, while those at A = l were obtained from 
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the fit to Lanczos results carried out in Ref. [|L3|. We observe that for A<0.5 the agreement 



with exact results is relatively good even for the states at the top of the one-magnon band. 

We now proceed to discuss the two dynamic spin responses. In Fig. f| we display our 
results for S (q, uj) and for S n (q, uj) as a function of the excitation energy uj for various values 
of the momentum transfer q. For clarity, all of the So responses have been multiplied by a 
factor of ten, while the S w responses include an artificial smearing of the delta function: 

5(uj-uj )^- ^ - ; r? = 0.02. (52) 

As suggested in our earlier discussion the S response is indeed small and — as it probes the 
two-magnon structure of the system — it is located at high-excitation energy and its strength 
is strongly fragmented. In contrast, the S n response is concentrated in a single excitation 
and, at least for q ~ tt, it is located at low-excitation energy (u ^ A). The fact that all 
the S n strength is concentrated in one fragment suggests that, once the dispersion band 
is computed [see. Fig. |3|], all that remains to fully characterize the S n (q, uj) response is to 
determine the static structure factor S n (q). Such a plot is displayed in Fig. |5] where So(q) 
also appears. The figure clearly shows that 5*0 (?) is not only appreciably smaller than S n (q) 
but it is also considerably smoother. This suggests that the large Fourier components of 
the sum of the two structure factors — and thus the long-range behavior of the spin-spin 
correlation function — will, indeed, be dominated by S n (q). Thus, Eq. (|40| ) can be used to 
estimate the spin-spin correlation length: 

£ « —r 1 1 — ► = 2.323 . (53) 



In 



d/(l- vT^d 5 ) 



d=0.914 



Note that a value of £ = 2.211 for the correlation length, extracted from an exponential 
fit to the full numerical evaluation of the spin-spin correlation function, compares well with 
the simple analytic estimate given above — but differs substantially from the exact value of 
£ = 3.19§. 



B. Extension to Larger Bases 

In this section we compare mean-field results obtained on a variety of angular-momentum 
coupled bases. We have alluded earlier to the fact that for gapped systems with a relatively 
small correlation length — such as the two-leg ladder materials studied here — incorporating 
as much as possible of the short-distance structure of the system into the definition of the 
basis is clearly advantageous [|,14]- The three angular-momentum coupled bases considered 



here are defined as the eigenstates of the 2 x L Heisenberg ladder, with L — 1, 2, and 4; 
for L = 1 and 2 these are the rung and plaquette basis, respectively. For simplicity — 
and because of their intrinsic interest — we limit ourselves to the physically relevant case of 
isotropic (J± = Jn = 1) ladders. 

In Table |l| we list results for the ground-state energy-per-site and for the single- 
triplet gap in various approximations; the next to last column contains the exact — basis- 
independent — results obtained from the DMRG calculation of Ref. || , while the last column 
lists the spin-singlet condensate fraction. This last column makes evident that the "purity" 
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of the ground-state increases as one incorporates more and more of the short-range structure 
of the system. However, as severe truncations must be made, the overall picture does not 
necessarily improve with the size of the basis. The second column contains the zeroth-order 
results. These results are obtained from the mere definition of the basis and are independent 
of any approximation. We observe that for these gapped systems incorporating as much as 
possible of the short-range dynamics into the definition of the basis is, indeed, advantageous. 
For example, for the 2x4 basis the zeroth-order result for the energy is already within seven 
percent of the exact answer. For the gap — a much more sensitive quantity — the result still 
differs from the exact answer by 50 percent. Yet, this is significantly better than the factor- 
of-two discrepancy obtained in the rung and plaquette bases. The third column lists the 
mean-field results. We observe that for the energy-per-site the reduction relative to the 
zeroth-order values is in good agreement with the second-order estimate from Eq. (|44a|) . 
That is, 

r -0.18750 for the 2 x 1 basis; 
[A£ /site] MF - { -0.04167 for the 2 x 2 basis; (54) 
I -0.01054 for the 2 x 4 basis. 

It is interesting to note that although the zeroth-order estimate for the energy-per-site varies 
substantially among the three bases (up to 30 percent) the mean-field approximation brings 
all of them to within one percent of each other — and to within approximately five percent of 
the exact answer. In the fourth column we have included the first-order correction coming 
from H 2 (see Sec.|ll E|): 



-0.04564 for the 2 x 1 basis; 
[A£ /site] H2 = { -0.00386 for the 2 x 2 basis; (55) 
-0.00050 for the 2 x 4 basis. 

In all cases the correction to the energy-per-site is small relative to the mean-field contribu- 
tion, and justifies — albeit a posteriori — treating H 2 perturbatively. 

The situation changes dramatically, however, upon examination of the singlet-triplet 
gap. In contrast to the energy-per-site which has a lowest-order contribution proportional 
to A 2 , the mean- field contribution to the spin gap starts with the term — 2A|a|. Although, 
for A = 1, the correction from the higher-order terms is not negligible, this linear term 
dominates. Note that while the mean-field correction to the spin gap is close to 40 percent 
in the 2x4 basis, the corresponding correction to the energy-per-site is a meager 2 percent. 
Moreover, it is now no longer justifiable — specifically in the case of the rung basis — to treat 
H 2 perturbatively. Indeed, the almost-perfect agreement between the rung-basis mean-field 
and the exact answer is lost. While, on its own, this result does not invalidate the mean- field 
approach, a careful re-assessment of the role of H 2 is clearly necessary. 

The impact of H 2 on the gap, however, seems to be more moderate for the larger bases. 
For example, there is only a ten and three percent correction to the gap in the plaquette and 
2x4 basis, respectively. This seems to indicate that the whole approach, namely, mean-field 
plus H 2 , is under control for these larger bases. Yet the agreement with the exact result in 
the particular case of the 2x4 basis is poor (of the order of 15 percent). We attribute this 
behavior to the severe truncation of the basis; recall that in the 2x4 basis only four out of 
a total of 256 states have been retained. 
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Overall, it seems that the plaquette basis offers the best compromise. While there are 
non trivial perturbative corrections to the gap from H 2 , these are moderate — unlike the 
corrections computed in the rung basis. Moreover, while we do truncate the basis, the 
truncation errors are not as severe as in the 2x4 basis. This is reflected in the "bracketing" 
of the exact result for the gap between the mean-field and the mean-field plus H2 estimates. 
Further, because of the "small" expansion parameter in the plaquette basis [(Aa/Ao) 2 = 
A 2 /9] most of the mean-field behavior can be reproduced with terms up to order A 2 in the 
perturbative expansion. This is illustrated in Fig. ^| where we have plotted the ground-state 
energy-per-site and the spin gap as a function of A. The filled circles, essentially an exact 
calculation, are from a previous perturbative estimate ||. Note that the A 7^ 1 case represents 
a pair of dimerized chains connected to each other by transverse bonds of strength equal to 
one of the two intra-chain couplings. To the best of our knowledge, a physical realization 
of this "dimerized ladder" is yet to be found. Yet, as it is well known, the isotropic limit of 
these ladders has been the object of intense scrutiny [TTI . 



We conclude this section with a discussion of plaquette-basis results for the dynamic spin 
response S n (q, u); recall that in a mean- field approximation with the low-energy (truncated) 
plaquette basis the So(q,u) response is identically zero. As the strength of the S n response 
is concentrated in a single excitation, it is sufficient to present results for the one-magnon 
dispersion band and for the static structure factor in order to fully characterize its behavior. 
Thus, in Fig. [7] we display the dispersion relation for the one-magnon band at the isotropic 
value. The solid and dashed lines (practically indistinguishable in the figure) represent the 
mean-field and the mean-field plus H2 results, respectively. The filled circles were computed 
from a a fit to Lanczos data JIB] . We observe that the plaquette results describe accurately 
the high- momentum/low-energy part of the band, but deteriorate rapidly as one moves to 
the low-momentum/high-energy region. This behavior is to be expected, as the discarded 
(high-energy) states in the plaquette basis account for most of the high-energy properties of 
the systems. In the near future we are planning to incorporate additional high-energy states 
into the mean-field model in order to improve the description of the high-energy properties 
of the system. Fortunately, as is evident in Fig. [| [see also Fig. ||, the discarded high-energy 
states are of little consequence to the dynamic spin response, which has most of its strength 
concentrated at low excitation energy. Moreover, this is the region that controls the long- 
range behavior of the spin-spin correlation function, which is shown in the inset on Fig. |8|. 
According to Eq. ([51]), the plaquette-basis spin-spin correlation function per site is then 
given by 

<*„| S ,(0) S , W |* > = t f ? : ■ («) 

l + dcos(2g) 



6 Jo 2tt 

A steepest-descent evaluation of the integral on the rhs of this equation yields, in analogy 
with Eq. ( |40|) for the rung basis, 

(%\S z (0)S z (r)\%) ~ e"^ ; £ = 2/ In [-d/(l - VT^)} (57) 

where we must recall that, in the plaquette basis, d<0, while d is positive in the rung basis. 
A fit to the full result for the correlation function (see Fig. (§)) gives a value of £ = 3.098 for 
the correlation length, which is within seven percent of the approximate analytic value of 
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These results compare very well — at the three-percent level — with the exact value of £ = 3.19 
obtained from the DMRG analysis of Ref. ||. Our results for correlation lengths using both 
rung and plaquette bases are summarized in Table [TV. 



IV. CONCLUSIONS 

We have computed magnetic properties of Heisenberg spin- 1/2 ladders using an extension 
of the mean-field approximation of Gopalan, Rice, and Sigrist Q] which we have implemented 
using a variety of angular-momentum coupled bases. In the original version of the model — 
which employed the rung basis exclusively — the Heisenberg Hamiltonian was written in 
terms of spin-singlet and spin-triplet operators. Because of the rotational symmetry of the 
Hamiltonian, the full theory is specified in terms of only four parameters; these are the 
rung-singlet and rung-triplet energies, and two off-diagonal matrix elements. In the limit 
in which the ratio (A) of the coupling of the spins along the chains of the ladder relative to 
the coupling along the rungs is negligible, the ground-state wavefunction becomes a direct 
product of spin-singlets on every rung of the ladder. The first spin excitation is reached 
by breaking a spin-singlet along any of the rungs of the ladder. This weak-coupling limit 
of the theory is the motivation behind the mean- field approximation of Ref. JF|. Hence, 
in this approximation, one assumes that the spin-singlet bosons of zero linear momentum 
are condensed. This enables one to separate the Heisenberg Hamiltonian into a part that 
is quadratic in the spin-triplet operators t — the mean-field Hamiltonian — and a part {H 2 ) 
that remains quartic in t. The mean-field Hamiltonian is then brought to a diagonal form 
by means of a Bogoliubov transformation while the quartic part must be treated by some 
other means. 

As the mean-field Hamiltonian is determined by only a small number of parameters, its 
generalization to larger bases is straightforward. The physical motivation behind the change 
of basis is the small value (of the order of three lattice sites) of the spin-spin correlation 
length. Thus, the main advantage of the larger bases is that they can incorporate exactly — by 
their mere definition — part of the important short-range dynamics of these gapped systems. 

Our first calculations were carried out in the rung basis, as in Ref. . We implemented 
the mean-field procedure through a Bogoliubov transformation, which although standard, 
was effected here via a diagonalization of a suitable RPA-like matrix. One of the virtues of 
such an approach is that the generalization to larger vector spaces becomes simple. This 
is of special interest to us, as our future mean-field work will incorporate additional high- 
energy states of the plaquette basis. Our rung-basis results differ from those presented by the 
authors of Ref. [0] due to what we believe was an oversight on their part. Our results indicate 
that, contrary to their first reports, their mean-field approach works well at a quantitative 
level. Indeed, we find that the results for the energy-per-site and for the singlet-triplet gap 
agree with exact calculations to within six percent all through the < A < 1 range (note, 
however, that the value for the spin-spin correlation length at the isotropic point was heavily 
underestimated: £ = 2.211 < 3.19). In particular, the spin gap differs from the exact answer 
by less than one percent at the isotropic point. Unfortunately, this nice agreement was lost 
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once the contribution from H 2 was estimated perturbatively. Incorporating the effects from 
H 2 in a reliable manner should constitute an important goal for future work. 

For the larger 2x2 (plaquette) and 2x4 bases considered here an additional approxi- 
mation was required. We assumed that the important physics of the spin-ladder materials 
is dominated by the lowest singlet and triplet states in the corresponding bases; note that it 
is only in the rung basis that the singlet and triplet states are unique. We found that while 
rung- and plaquette-basis results are of similar quality for most ladder observables, only the 
latter ones were stable upon including the perturbative contribution from H 2 . This is a 
direct consequence of the size of the elementary block defining the basis. As the larger bases 
incorporate more of the short-range dynamics of the system, the quartic piece of the Hamil- 
tonian (H2) decreases in importance as the size of the elementary block increases. Since 
our general approach is easily extended to any angular-momentum coupled basis, mean-field 
studies were also carried out utilizing a 2 x 4 basis. The results obtained with this basis did 
not improve — and in some cases worsen — the agreement relative to the plaquette basis. We 
have traced this behavior to the severe truncations made in the 2x4 basis; only four out of a 
total of 256 states were retained. Thus, we conclude that the plaquette basis offers the best 
compromise between building enough short-range physics into the definition of the basis, 
while avoiding too severe a truncation. Indeed, plaquette-basis mean-field results for the 
energy-per-site, the spin gap, and the spin-spin correlation function are within five percent 
of the quoted DMRG values 0. 

In the future, we plan to complete a number of extensions of the mean-field treatment. 
For example, we have truncated the spectrum of spin triplet states in the plaquette ba- 
sis retaining only those of lowest energy. Although we have observed that this is a good 
approximation for describing the low-energy dynamics, it is not adequate to describe the 
high-energy properties. This requires an extension of the usual Bogoliubov technique for 
bringing quadratic Hamiltonians involving a single species of field operator into diagonal 
form. We have already made considerable progress in this direction and we are planning 
to publish our results in the near future. Alternatively, one could still retain the truncated 
basis but account for the truncation errors through a suitable renormalization of the ma- 
trix elements. Indeed, one such technique — CORE, or the contractor renormalization group 



method |15] — has been used successfully by us in computing static and dynamic properties 
of up to 2 x 16 ladders ||. Such an approach, to be used in the near future, could also be 
valuable in the context of the mean-field theory. We are also planning to develop a version 
of the mean-field theory for ladders which treats the quartic terms of the Hamiltonian self- 
consistently, rather than just perturbatively as it was done here. Finally — and of special 
relevance to the high-T c materials — we are planning to use the mean-field theory in the 



various angular-momentum coupled bases to study the physics of lightly-doped ladders [16 



APPENDIX: BOGOLIUBOV TRANSFORMATION 

In this appendix we carry out the Bogoliubov transformation via a diagonalization of 
a random-phase-approximation (RPA) matrix. In order to do so we write the — operator 
part — of the mean-field Hamiltonian of Eq. (F2I|) in matrix form. That is, 
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K UF (s 2 ,fi) = l -Y. T m {k)mT m {k) . (Al) 

km 

Although we illustrate the method using the 2x2 matrix of Eq. (|22[) the method can be 
easily generalized to matrices of arbitrary dimension. 

Given that the matrix Q(k) is symmetric, it can be diagonalized via a similarity trans- 
formation. The problem with such an approach, however, is that the resulting orthogonal 
matrix generates a new set of creation and annihilation operators that are not consistent 
with each other. Instead, what is required is to write Eq. ( |A~1| ) in terms of a matrix having 
an RPA-like structure, i.e., 

K MF (s 2 ,fi) = ^Y, T i(k)gn RPA (k)T m (k) , (A2) 

km 

where we have introduced the "metric" g = diag{l, —1} and have defined Qrpa = g^- We 
now list, without proof, some useful properties of RPA matrices: 

1. Eigenvalues of Qrj>a come in pairs; if ui is an eigenvalue of Qrpa so is — to. 

2. Eigenvectors of ^rpa corresponding to different eigenvalues are orthonormal with re- 
spect to the metric g: (u'\g\cj) = sgn(u;) 

3. There exists a "similarity" transformation S — with its columns being the eigenvectors 
of fi RP A — that brings Qrpa into a diagonal form: 5 -1 f2 RP AiS = diagju;, — u} , where 
S- 1 =gS T g. 

Using these properties, and defining a new set of creation and annihilation operators through 
T m = S~ l T m , the mean-field Hamiltonian becomes explicitly diagonal. That is, 

M^ri^ErlWft' w ^)r m (A;) = |E^)+E«(*)7!L(*)7m(A:). (A3) 

km \ J ' k k m 
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FIG. 1. Ground-state energy-per-site (in units of J±) as a function of A. The filled circles are 
"exact" results obtained from a very high-order perturbative expansion ||], while GRS are the 
results of Gopalan, Rice, and Sigrist [J7|. 
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FIG. 2. Singlet-triplet gap (in units of J±) as a function of A. The filled circles are "exact" 
results obtained from a very high-order perturbative expansion Q, while GRS are the results of 
Gopalan, Rice, and Sigrist Q. 
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FIG. 3. Spin-triplet dispersion band relative to the band minimum (in units of J±) as a function 
of momentum. The filled symbols are exact results H|i~3f (see text); the dotted lines joining them 
were obtained from a cubic-spline fit. 
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FIG. 4. Dynamic spin responses So and S n in the rung basis as a function of lo for various 
values of q. The responses were computed at the isotropic value of A = 1. All Sq responses have 
been multiplied by a factor of 10, while the S n responses include a smearing factor of 0.02. 
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FIG. 5. Static structure factors So and in the rung basis as a function of q for A = 1. The 
inset shows a logarithmic plot of the spin-spin correlation as a function of r. A value of £ = 2.211 
is obtained for the spin-spin correlation length from the slope of the straight-line fit. 
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FIG. 6. Energy-per-site and singlet-triplet gap as a function of A in the plaquette basis. Quan- 
tities plotted are ratios to their exact values at the isotropic (A = 1) point. The filled circles are 
from high-order perturbative calculations M. 
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FIG. 8. Static structure factors S n in the plaquette basis as a function of q at the isotropic 
point; the structure factor So vanishes in the (low-energy) plaquette basis. The inset shows a 
logarithmic plot of its Fourier transform; the indicated straight-line fit yields a value of £ = 3.098 
for the spin-spin correlation length. 



27 



TABLES 



TABLE I. Unperturbed singlet/triplet energies and the two independent matrix elements in 
various angular-momentum-coupled bases. 

Basis eo ei a = (011||V(1,2)||101) (3 = (110||V(1,2)||110) 

2 x 1 -3/4 +1/4 +1/2 -1 

2x2 -2 -1 -1/3 -1/4 

2 x 4 -4.29307 -3.52286 -0.20805 -0.08625 



TABLE II. Results for the ground-state energy-per-site and for the singlet-triplet gap (dis- 
played inside brackets) for a few values of A in the Rung basis. 

A Ref. (?J Present work Exact 

025 -0.37908 [0.72846] -0.38683 [0.76699] -0.38804 [0.78438] 

0.50 -0.39392 [0.40411] -0.42237 [0.59685] -0.42911 [0.64065] 

1.00 -0.47531 [0.10291] -0.54285 [0.50133] -0.57804 [0.50400] 



TABLE III. Results for the ground-state energy-per-site and for the singlet-triplet gap (dis- 
played inside brackets) for various angular-momentum coupled bases. Listed are zeroth-order, 
mean- field, mean-field plus H%, and exact results (see text for details). The exact values were 
extracted from DMRG calculations. The last column lists the spin-singlet condensate fraction. 

Basis Zeroth-order MFT MFT+H 2 Exact s 2 

2 x 1 -0.37500 [1.00000] -0.54285 [0.50133] -0.58849 [0.13204] -0.57804 [0.50400] 0.78129 

2 x 2 -0.50000 [1.00000] -0.54112 [0.53689] -0.54498 [0.48154] -0.57804 [0.50400] 0.85826 

2 x 4 -0.53663 [0.77021] -0.54724 [0.44566] -0.54774 [0.43216] -0.57804 [0.50400] 0.89628 



TABLE IV. Results for spin spin-spin correlation lengths for the rung and plaquette bases. 
The "Fit" results are obtained by straight-line fits to semi-log plots of the correlation function as 
determined by numerical evaluation of the integral appearing in Eqn.(^7|). The "Analytic" results 

for rung and plaquette bases, respectively. The "Exact" 

H. 



are obtained using Eqns.(53) and 
value is from the DMRG calculation of Ref. 



Basis 


Fit 


Analytic 


Exact 


2 x 1 


2.211 


2.323 


3.19 


2x2 


3.098 


3.327 


3.19 



28 



